Accessing Data Using MODIStsp
Aditya Ranganath
3/20/2022
library(MODIStsp)
library(sf)
library(tidyverse)
library(tmap)
library(tidycensus)
library(raster)
library(viridis)In order to use the MODIStsp package, you must register on the EARTHDATA page, at this link.
1 MODIS Land Cover Data
We will extract the MODIS land cover raster (MCD12C1) and make a simple map of land cover in the state of Colorado. A description of the data can be found here.
1.1 Extract and download Colorado polygon
First, let’s extract a polygon of Colorado’s state borders using tidycensus.
co_polygon<-get_decennial(geography = "state",
variables = "P001001",
year = 2010,
geometry=TRUE) %>%
filter(NAME=="Colorado") %>%
st_transform(4326)In case we want to overlay county-level boundaries against the land cover raster, we can also extract a spatial dataset of Colorado counties:
co_counties_polygon<-get_decennial(geography = "county",
state="CO",
variables = "P001001",
year = 2010,
geometry=TRUE) %>%
st_transform(4326)After extracting this data from tidycensus, we’ll need to extract its bounding box, or alternatively, download the data to our local machine; this will allow us to specify the desired extent of the raster when we extract it using the MODIStsp function. In this case, we’ll go with the first option.
First, we’ll define a path and name for the spatial dataset that is to be downloaded:
filepath<-"land_cover/colorado.shp"Then, we’ll download the co_polygon spatial data object as a shapefile using the st_write function:
st_write(co_polygon, paste0(filepath))The Colorado state polygon is stored in a sub-directory (named “land_cover”) within our working directory, and is named “colorado.shp”.
1.2 Extract land cover raster from MODIS
Now, we’re ready to
MODIStsp_get_prodlayers("MCD12Q1")## $prodname
## [1] "LandCover_Type_Yearly_500m (MCD12Q1)"
##
## $bandnames
## [1] "LC1" "LC2" "LC3"
## [4] "LC4" "LC5" "LC_Prop1"
## [7] "LC_Prop2" "LC_Prop3" "LC_Prop1_Assessment"
## [10] "LC_Prop2_Assessment" "LC_Prop3_Assessment" "LC_QC"
## [13] "LC_LW"
##
## $bandfullnames
## [1] "Land Cover Type 1 (IGBP)*"
## [2] "Land Cover Type 2 (UMD)*"
## [3] "Land Cover Type 3 (LAI/fPAR)*"
## [4] "Land Cover Type 4 (NPP/BGC)*"
## [5] "Land Cover Type 5: Annual Plant Functional Types classification"
## [6] "FAO-Land Cover Classification System 1 (LCCS1) land cover layer"
## [7] "FAO-LCCS2 land use layer"
## [8] "FAO-LCCS3 surface hydrology layer"
## [9] "LCCS1 land cover layer confidence"
## [10] "LCCS2 land use layer confidence"
## [11] "LCCS3 surface hydrology layer confidence"
## [12] "Land Cover QC"
## [13] "Binary land/water mask derived from MOD44W"
##
## $quality_bandnames
## NULL
##
## $quality_fullnames
## NULL
##
## $indexes_bandnames
## NULL
##
## $indexes_fullnames
## NULLMODIStsp(gui = FALSE,
out_folder = "land_cover",
out_folder_mod = "land_cover",
selprod = "LandCover_Type_Yearly_500m (MCD12Q1)",
bandsel = "LC1",
user = "YOUR_USERNAME" ,
password = "YOUR_PASSWORD",
start_date = "2019.01.01",
end_date = "2019.12.31",
output_proj = 4326,
verbose = FALSE,
spatmeth = "file",
spafile = filepath,
out_format = "GTiff")1.3 Read in land use raster
CO_raster<-raster("land_cover/colorado/LandCover_Type_Yearly_500m_v6/LC1/MCD12Q1_LC1_2019_001.tif")1.4 Visualize raster with tmap
raster_label<-c("Evergreen needleleaf forests",
"Evergreen broadleaf forests",
"Deciduous needleleaf forests",
"Deciduous broadleaf forests",
"Mixed forests",
"Closed shrublands",
"Open shrublands",
"Woody savannas",
"Savannas",
"Grasslands",
"Permanent wetlands",
"Croplands",
"Urban and built-up lands",
"Cropland/natural vegetation mosaics",
"Snow and ice",
"Barren",
"Water bodies")CO_landcover_tmap<-tm_shape(CO_raster)+
tm_raster(breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),
labels=raster_label,
palette="Set1",
title="Colorado Land Cover, 2019")+
tm_legend(outside=T)tmap_mode("plot")## tmap mode set to plottingCO_landcover_tmap## Some legend labels were too wide. These labels have been resized to 0.58. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tmap_mode("view")## tmap mode set to interactive viewingCO_landcover_tmaptmap_mode("plot")## tmap mode set to plottingCO_landcover_tmap+
tm_shape(co_counties_polygon)+
tm_polygons(alpha=0)## Some legend labels were too wide. These labels have been resized to 0.58. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tmap_mode("view")## tmap mode set to interactive viewingCO_landcover_tmap+
tm_shape(co_counties_polygon)+
tm_polygons(alpha=0)1.5 Visualize raster with ggplot
CO_raster_df <- as.data.frame(CO_raster, xy = TRUE) %>%
mutate(MCD12Q1_LC1_2019_001 = as.factor(MCD12Q1_LC1_2019_001))levels(CO_raster_df$MCD12Q1_LC1_2019_001)## [1] "1" "4" "5" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17"CO_raster_df$MCD12Q1_LC1_2019_001<-factor(CO_raster_df$MCD12Q1_LC1_2019_001, levels=c("1", "2",
"3", "4", "5",
"6", "7", "8",
"9", "10", "11",
"12", "13", "14",
"15", "16", "17"))levels(CO_raster_df$MCD12Q1_LC1_2019_001) <- raster_labellevels(CO_raster_df$MCD12Q1_LC1_2019_001)## [1] "Evergreen needleleaf forests" "Evergreen broadleaf forests"
## [3] "Deciduous needleleaf forests" "Deciduous broadleaf forests"
## [5] "Mixed forests" "Closed shrublands"
## [7] "Open shrublands" "Woody savannas"
## [9] "Savannas" "Grasslands"
## [11] "Permanent wetlands" "Croplands"
## [13] "Urban and built-up lands" "Cropland/natural vegetation mosaics"
## [15] "Snow and ice" "Barren"
## [17] "Water bodies"ggplot() +
geom_raster(data = CO_raster_df,
aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
geom_sf(data = co_polygon, inherit.aes = FALSE, fill = NA) +
scale_fill_viridis(name = "Land Cover Type", discrete = TRUE) +
labs(title = "Land Cover classification in Colorado",
subtitle = "01-01-2019 - 31-12-2019",
x = "Longitude",
y = "Latitude") +
theme_minimal()
ggplot() +
geom_raster(data = CO_raster_df,
aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
geom_sf(data=co_counties_polygon, inherit.aes=F, fill=NA)+
scale_fill_viridis_d(name = "Land Cover Type", option="C", direction=-1 ) +
labs(title = "Land Cover classification in Colorado",
subtitle = "01-01-2019 - 31-12-2019",
x = "Longitude",
y = "Latitude") +
theme_minimal()
CO_landcover_ggplot<-
ggplot() +
geom_raster(data = CO_raster_df,
aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
geom_sf(data=co_counties_polygon, inherit.aes=F, fill=NA)+
labs(title = "Land Cover classification in Colorado",
subtitle = "01-01-2019 - 31-12-2019",
x = "Longitude",
y = "Latitude") +
theme_minimal()colors<-c("lightpink1", "lightgrey", "magenta", "lightcyan",
"yellow", "palegreen", "seashell", "plum",
"skyblue", "salmon2", "khaki", "tomato",
"gray0", "springgreen", "purple", "orange",
"mediumblue", "slategray")
CO_landcover_ggplot +
scale_fill_manual(values=colors, drop=F)
2 Extracting and visualizing an NDVI raster for Colorado
MODIStsp_get_prodlayers("M*D13A2")## $prodname
## [1] "Vegetation_Indexes_16Days_1Km (M*D13A2)"
##
## $bandnames
## [1] "NDVI" "EVI" "VI_QA" "b1_Red" "b2_NIR" "b3_Blue"
## [7] "b7_SWIR" "View_Zen" "Sun_Zen" "Rel_Az" "DOY" "Rely"
##
## $bandfullnames
## [1] "16 day NDVI average" "16 day EVI average"
## [3] "VI quality indicators" "Surface Reflectance Band 1"
## [5] "Surface Reflectance Band 2" "Surface Reflectance Band 3"
## [7] "Surface Reflectance Band 7" "View zenith angle of VI pixel"
## [9] "Sun zenith angle of VI pixel" "Relative azimuth angle of VI pixel"
## [11] "Day of year of VI pixel" "Quality reliability of VI pixel"
##
## $quality_bandnames
## [1] "QA_qual" "QA_usef" "QA_aer" "QA_adj_cld" "QA_BRDF"
## [6] "QA_mix_cld" "QA_land_wat" "QA_snow_ice" "QA_shd"
##
## $quality_fullnames
## [1] "VI Quality"
## [2] "VI usefulness"
## [3] "Aerosol quantity"
## [4] "Adjacent cloud detected"
## [5] "Atmosphere BRDF correction performed"
## [6] "Mixed Clouds"
## [7] "Land/Water Flag"
## [8] "Possible snow/ice"
## [9] "Possible shadow"
##
## $indexes_bandnames
## [1] "SR" "NDFI" "NDII7" "SAVI"
##
## $indexes_fullnames
## [1] "Simple Ratio (NIR/Red)"
## [2] "Flood Index (Red-SWIR2)/(Red+SWIR2)"
## [3] "NDII7 (NIR-SWIR2)/(NIR+SWIR2)"
## [4] "SAVI (NIR-Red)/(NIR+Red+0.5)*(1+0.5)"MODIStsp(
gui = FALSE,
out_folder = "VegetationData",
out_folder_mod = "VegetationData",
selprod = "Vegetation_Indexes_16Days_1Km (M*D13A2)",
bandsel = "NDVI",
user = "YOUR_USERNAME",
password = "YOUR_PASSWORD",
start_date = "2020.07.03",
end_date = "2020.07.03",
output_proj = 4326,
verbose = FALSE,
spatmeth = "file",
spafile = "VegetationData/colorado.shp",
out_format = "GTiff"
)2.1 Load NDVI raster
CO_ndvi<-raster("VegetationData/colorado/VI_16Days_1Km_v6/NDVI/MYD13A2_NDVI_2020_185.tif")CO_ndvi## class : RasterLayer
## dimensions : 269, 471, 126699 (nrow, ncol, ncell)
## resolution : 0.01490176, 0.01491085 (x, y)
## extent : -109.0603, -102.0415, 36.99243, 41.00344 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /Users/adra7980/Documents/git_repositories/ARSC5040_GIS/general_resources/VegetationData/colorado/VI_16Days_1Km_v6/NDVI/MYD13A2_NDVI_2020_185.tif
## names : MYD13A2_NDVI_2020_185
## values : -32768, 32767 (min, max)See here for scale factor.
CO_ndvi_adjusted<-CO_ndvi*0.00012.2 Make ndvi map using tmap
tmap_mode("plot")## tmap mode set to plottingCO_ndvi_map<-tm_shape(CO_ndvi_adjusted)+
tm_raster(palette = "YlGn", style="cont", breaks=c(-0.2, 0,0.2,0.4,0.6,0.8, 1))+
tm_legend(outside=T)
CO_ndvi_map## Warning: Breaks contains positive and negative values. Better is to use
## diverging scale instead, or set auto.palette.mapping to FALSE.
tmap_mode("view")## tmap mode set to interactive viewingCO_ndvi_map## Warning: Breaks contains positive and negative values. Better is to use
## diverging scale instead, or set auto.palette.mapping to FALSE.2.3 Make ndvi map using ggplot
NDVI_co_df <- as.data.frame(CO_ndvi_adjusted, xy = TRUE, na.rm = TRUE)
# Visualising using ggplot2
ggplot() +
geom_raster(
data = NDVI_co_df,
aes(x = x, y = y, fill = MYD13A2_NDVI_2020_185)
) +
geom_sf(data = co_counties_polygon, inherit.aes = FALSE, fill = NA) +
scale_fill_viridis(name = "NDVI") +
labs(
title = "NDVI (Normalized Difference Vegetation Index) in Colorado",
subtitle = "03-07-2020",
x = "Longitude",
y = "Latitude"
) +
theme_minimal()